% function test_biharmonic_h
clear all;
global V posV T posT E posE TE ET;

caseNum = 3; disp_d = 7; tri_temp = make_mesh(disp_d); 

[p,e,t] = initmesh('lshapeg','Hmax',0.4,'Hgrad',1.99);
V = [p(1,:)' -p(2,:)']; T = t(1:3,:)';
posV = size(V,1); posT = size(T,1); % init posV and posT here

% reserve more spaces
V = [V;zeros(10000,2)]; T = [T;zeros(10000,3)];
T_more = zeros(size(T,1),3); % additional information for elements
d = 5*ones(size(T,1),1);

max_err = 1e-5;   % 1e-3 1e-4 1e-5
refine_L_shape(10); mkcc; % 8 9 10
posE = 0; areas = tdata; % init and reserve space for E,TE,ET

% prepare the global bending matrix for saving time
% load com15 cr desc asce Int_UU_cell Int_UV_cell Int_UW_cell Int_VV_cell Int_VW_cell Int_WW_cell Int_Mass MatVDM;
d_max = max(d);  [cr,desc,asce] = get_auxillary_mat(d_max);
% value for assemble
Int_UU_cell = cell(d_max,1);
Int_UV_cell = cell(d_max,1);
Int_UW_cell = cell(d_max,1);
Int_VV_cell = cell(d_max,1);
Int_VW_cell = cell(d_max,1);
Int_WW_cell = cell(d_max,1);
Int_Mass = cell(d_max,1);
MatVDM = cell(d_max,1);
for degree = 5:d_max
    Int_Mass{degree} = build(degree);
    MatVDM{degree} = vdm22(degree,disp_d);
    Mat = build(degree-2);
    Id = diag(ones((degree+1)*(degree+2)/2,1));
    Du = degree*de_cast_step(Id,degree,1,0,-1,desc); % the direction derivate on v1-v3
    Duu = (degree-1)*de_cast_step(Du,degree-1,1,0,-1,desc); % the direction derivate on v1-v3
    Duv = (degree-1)*de_cast_step(Du,degree-1,0,1,-1,desc); % blending partial derivatives
    Dv = degree*de_cast_step(Id,degree,0,1,-1,desc);  % the direction derivate on v2-v3
    Dvv = (degree-1)*de_cast_step(Dv,degree-1,0,1,-1,desc); % the direction derivate on v2-v3
    Int_UU_cell{degree} = Duu'*Mat*Duu;
    Int_UV_cell{degree} = Duu'*Mat*Dvv;
    Int_UW_cell{degree} = Duu'*Mat*Duv;
    Int_VV_cell{degree} = Dvv'*Mat*Dvv;
    Int_VW_cell{degree} = Dvv'*Mat*Duv;
    Int_WW_cell{degree} = Duv'*Mat*Duv;
end
% value for posteriori estimate:
[qw,qp] = quad_tri(13);
qw = qw'; qx = qp(1,:)'; qy = qp(2,:)';
phi4 = cell(d_max,5);
for degree = 5:d_max
    Id = diag(ones((degree+1)*(degree+2)/2,1));
    
    Du = degree*de_cast_step(Id,degree,1,0,-1,desc); % the direction derivate on v1-v3
    Dv = degree*de_cast_step(Id,degree,0,1,-1,desc); % the direction derivate on v1-v3
    
    Duu = (degree-1)*de_cast_step(Du,degree-1,1,0,-1,desc); 
    Duv = (degree-1)*de_cast_step(Du,degree-1,0,1,-1,desc); 
    Dvv = (degree-1)*de_cast_step(Dv,degree-1,0,1,-1,desc); 
    
    Duuu = (degree-2)*de_cast_step(Duu,degree-2,1,0,-1,desc);
    Duuv = (degree-2)*de_cast_step(Duu,degree-2,0,1,-1,desc);  
    Duvv = (degree-2)*de_cast_step(Duv,degree-2,0,1,-1,desc); 
    Dvvv = (degree-2)*de_cast_step(Dvv,degree-2,0,1,-1,desc); 
    
    Duuuu = (degree-3)*de_cast_step(Duuu,degree-3,1,0,-1,desc); 
    Duuuv = (degree-3)*de_cast_step(Duuu,degree-3,0,1,-1,desc); 
    Duuvv = (degree-3)*de_cast_step(Duuv,degree-3,0,1,-1,desc); 
    Duvvv = (degree-3)*de_cast_step(Duvv,degree-3,0,1,-1,desc); 
    Dvvvv = (degree-3)*de_cast_step(Dvvv,degree-3,0,1,-1,desc); 
    % evaluate the values on quadrature points of all basis' 4th derivatives:
    Mat = vdm23(degree-4,qx,qy,1-qx-qy);
    phi4{degree,1} = Mat*Duuuu;
    phi4{degree,2} = Mat*Duuuv;
    phi4{degree,3} = Mat*Duuvv;
    phi4{degree,4} = Mat*Duvvv;
    phi4{degree,5} = Mat*Dvvvv;
end   

tic; count = 1; flag = 1;  % control the computing
while flag
    flag = 0; % solved again , then keep the flag 0    
    [dofs,n_dof] = sort_dof_dis2(T,posT,d); % init and reserve spaces for dofs and n_dofs                
    % the system equations 
    fprintf('\n\nThe %d step: Assembling systems...%d element \n',count,posT);
    K = bending_bb_2(dofs,V,T,posT,d,areas,Int_UU_cell,Int_UV_cell,Int_UW_cell,Int_VV_cell,Int_VW_cell,Int_WW_cell);
    M = mass_bb_2(dofs,V,T,posT,d,areas,Int_Mass);
    F = bnet_2(dofs,V,T,posT,d,'biharmonic_f',caseNum);     
    b1 = M*F;
    [B,G] = biharmonic_bc_xy(dofs,d,'biharmonic_g','biharmonic_gx','biharmonic_gy',cr,caseNum);
    
    % the smoothness conditions
    H = cr_matrix(dofs,V,T,E,posE,TE,ET,desc,asce,cr,d,1);
   
    % solute the linear algebra system
    k = length(G);    p = size(H,1);    n = size(K,1);    
    fprintf('Solving system ... matrix size : %d x %d \n',n,n);
    c = lagrange22(K,b1,[B;H],[G;zeros(p,1)]);
    
    fprintf('\nComputing the L_inf error...'); % adaptive procedure
%     [eta_L2,err] = comput_err_easy(dofs,V,T,posT,d,c,'biharmonic_u',caseNum);
    err = posteriori_biharmonic(dofs,V,T,posT,d,c,qw,qx,qy,phi4);
    err = smooth_error(V,posV,T,posT,err);
    error = sum(err);
    fprintf('\nthe total error is %e. matrix size is: %d\n',error,size(K,1));
    if error >= max_err
        flag = 1; % need to comput once more
        fprintf('\nRefineing the mesh...');
        refine_T = mesh_refine_by_fraction(err,0.8,0);
        refine_E = mark_bisect(V,T,TE,posE,posT,refine_T);
        areas = refine_bisect(refine_E);
        count = count +1;
    end
end
trisurf(Tris,Points(:,1),Points(:,2),Z);shading interp; view(-10,30);
figure('Position',[1 250 800 300]);
subplot(1,3,1);plot_t(V,T(1:posT,:));
axis([min(V(1:posV,1)) max(V(1:posV,1)) min(V(1:posV,2)) max(V(1:posV,2))]);
subplot(1,3,2);plot_t(V,T(1:posT,:));
axis([min(V(1:posV,1)) max(V(1:posV,1)) min(V(1:posV,2)) max(V(1:posV,2))]);
subplot(1,3,3);plot_t(V,T(1:posT,:));
axis([min(V(1:posV,1)) max(V(1:posV,1)) min(V(1:posV,2)) max(V(1:posV,2))]);
toc;

figure;
[Z,Tris,Points] = heval(dofs,V,T,posT,c,d,tri_temp,disp_d);
exact = feval('biharmonic_u',Points(:,1),Points(:,2),caseNum);
trisurf(Tris,Points(:,1),Points(:,2),abs(Z-exact));
shading interp;view(50,30);axis([-1 1 -1 1 0 2e-8]);